module moments3jb
implicit none
contains
subroutine normjbmvu(MD,VD,matA,matB,b,eta,si,sigma,dmut,dvsig,pth)
use linear_operators
use matrixfuns
use testmomsrot
use testmoments2
use datahadler
integer, parameter:: n=3,p=6
integer, intent(in)::pth
double precision, intent(in)::b(n),eta,si,dmut(pth,n),dvsig(pth,n*n),sigma(n,n)
double precision, intent(out)::matA(pth,pth),matB(p,pth),MD(pth+p),VD(pth+p,pth+p)
double precision:: isigma(n,n),rsigma(n,n),irsigma(n,n),beta(n),dmach,tol&
&,matV(n,n),matL(n*(n+1)/2,n*n),matD(n*n,n*(n+1)/2),matDm(n*(n+1)/2,n*n),irsigmah(n,n)&
&,matmsk(n,p),matvsk(n*n,p),msk(p,1),beta1,nu,gam,mstar1(n,pth),mstar2(n*n,pth)&
&,m400,m004,m300,m120,me(n,1),dep(n,pth),drvsig(n*n,n*n)
integer i,j,irank

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

if (dot_product(beta(2:n),beta(2:n))<1.0d-10) then
	matV=eye(n)
else
!	Sigma=eye(n)
   matV(1,:)=(/-0.57735026918963d0, 0.81649658092773d0,0.0d0/)
   matv(2,:)=(/-0.57735026918963d0,-0.40824829046386d0,-0.70710678118655d0/)
   matv(3,:)=(/-0.57735026918963d0,-0.40824829046386d0, 0.70710678118655d0/)

!	matV(:,1)=beta/dsqrt(dot_product(beta,beta))
!	matV(1,2:3)=(/-0.50497514502648d0,-0.58533678384716d0/)
!	matV(2,2:3)=(/0.78236638609871d0,-0.0d0/)
!	matV(3,2:3)=(/-0.36455855607615d0,0.81079026232156d0/)

!	print*, "Please, supply orthogonal matrix!"
!	print*, "beta", beta
!	stop
	!matV(:,1)=beta/dsqrt(dot_product(beta,beta))
end if
irsigmah=matV.tx.irsigma

beta(1)=dsqrt(dot_product(beta,beta))
beta(2:3)=(/0.0d0,0.0d0/)

forall(i=1:pth)
MD(i)=0.0d0
end forall
call meanjbc(MD(pth+1:pth+p),beta,eta,si)


VD(1:pth,1:pth)=matmul(dmut,matmul(isigma,.t.dmut))&
&+.25d0*(((dvsig.x.kron(.t.irsigmah,.t.irsigmah))&
&.x.(emempkememp(beta,eta,si)-(vec(dble(eye(n))).xt.vec(dble(eye(n))))))&
&.x.(kron(irsigmah,irsigmah).xt.dvsig))&
&+.5d0*((dvsig.xt.kron(irsigmah,irsigmah)).x.(emkememp(beta,eta,si).x.(irsigmah.xt.dmut)))&
&+.5d0*(.t.(((dvsig.xt.kron(irsigmah,irsigmah)).x.(emkememp(beta,eta,si).x.(irsigmah.xt.dmut)))))

call varjbc(VD(pth+1:pth+p,pth+1:pth+p),beta,eta,si)



forall(i=1:n,j=1:p)
	matmsk(i,j)=0.0d0
end forall
beta1=dsqrt(dot_product(beta,beta))
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si

matmsk(1,1)=eubeta5(beta1,eta,si)
matmsk(1,2)=eub1hn(beta1,eta,si,2)*3.0d0
matmsk(1,3)=matmsk(1,2)
matmsk(1,4)=eubeta4(beta1,eta,si)
matmsk(2,5)=ehr(nu,gam,2)*3.0d0
matmsk(3,6)=matmsk(2,5)
forall(i=1:n*n,j=1:p)
	matvsk(i,j)=0.0d0
end forall
matvsk(1,1)=eubeta6(beta1,eta,si)-3.0d0
matvsk(5,1)=eub4hn(beta1,eta,si,1)-3.0d0
matvsk(9,1)=matvsk(5,1)


matvsk(1,2)=3.0d0*eub2hn(beta1,eta,si,2)-3.0d0
matvsk(5,2)=15.0d0*ehr(nu,gam,3)-3.0d0
matvsk(9,2)=3.0d0*ehr(nu,gam,3)-3.0d0

matvsk(1,3)=matvsk(1,2)
matvsk(5,3)=matvsk(9,2)
matvsk(9,3)=matvsk(5,2)

matvsk(1,4)=eubeta5(beta1,eta,si)
matvsk(5,4)=eub3hn(beta1,eta,si,1)
matvsk(9,4)=matvsk(5,4)

matvsk(2,5)=3.0d0*eub1hn(beta1,eta,si,2)
matvsk(4,5)=matvsk(2,5)

matvsk(3,6)=matvsk(2,5)
matvsk(7,6)=matvsk(2,5)

msk(:,1)=MD(pth+1:pth+p)


VD(1:pth,pth+1:pth+p)=matmul(dmut.xt.irsigmah,matmsk)&
&+.5d0*matmul(dvsig.xt.kron(irsigmah,irsigmah),matvsk-(vec(dble(eye(n))).xt.msk))

forall(i=1:pth,j=pth+1:pth+p)
VD(j,i)=VD(i,j)
end forall


matA=matmul(dmut,matmul(isigma,.t.dmut))&
&+.5d0*matmul(dvsig,matmul(kron(isigma,isigma),.t.dvsig))


forall(i=1:n*(n+1)/2,j=1:n*n)
	matL(i,j)=0.0d0
	matD(j,i)=0.0d0
end forall
matL(1:3,1:3)=dble(eye(3))
matL(4:5,5:6)=dble(eye(2))
matL(6,9)=1.0d0

matD(1:3,1:3)=dble(eye(3))
matD(4,2)=1.0d0
matD(7,3)=1.0d0
matD(5:6,4:5)=dble(eye(2))
matD(8,5)=1.0d0
matD(9,6)=1.0d0
matDm=matmul(.i.(matD.tx.matD),.t.matD)
mstar1=-(matmul(.t.matV,irsigma).xt.dmut)

drvsig=.5d0*((matL.tx.(.i.((matDm.x.kron(rsigma,dble(eye(n)))).xt.matL))).x.matDm)
mstar2=-kron(.t.matV,(.t.matV).x.irsigma)&
&.x.(drvsig.xt.dvsig)

call exportvec(vec2(drvsig),"drvsig.txt",2)


m400=eubeta4(beta1,eta,si)
m004=ehr(nu,gam,2)*3.0d0
m300=eubeta3(beta1,eta,si)
m120=eub1hn(beta1,eta,si,1)

!1
me(:,1)=(/m400,0.0d0,0.0d0/)
dep=(mstar1*m300+(kron(me,dble(eye(n))).tx.mstar2))
matB(1,:)=4.0d0*dep(1,:)
!2
me(:,1)=(/0.0d0,m004,0.0d0/)
dep=(kron(me,dble(eye(n))).tx.mstar2)
matB(2,:)=4.0d0*dep(2,:)
!3
me(:,1)=(/0.0d0,0.0d0,m004/)
dep=(kron(me,dble(eye(n))).tx.mstar2)
matB(3,:)=4.0d0*dep(3,:)
!4
me(:,1)=(/m300,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(4,:)=3.0d0*dep(1,:)
!5
me(:,1)=(/m120,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(5,:)=3.0d0*dep(2,:)
!6
me(:,1)=(/m120,0.0d0,0.0d0/)
dep=(mstar1+(kron(me,dble(eye(n))).tx.mstar2))
matB(6,:)=3.0d0*dep(3,:)

end subroutine normjbmvu


!!!
subroutine meanjbc(matmjb,beta,eta,si)
!Variance of the moment conditions of JB
use ghs
use testmomsrot
use linear_operators
double precision matmjb(6),beta(3),eta,si,beta1,nu,gam
integer i
beta1=dsqrt(dot_product(beta,beta))
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si

matmjb(1)=eubeta4(beta1,eta,si)-3.0d0
matmjb(2)=3.0d0*(ehr(nu,gam,2)-1.0d0)
matmjb(3)=matmjb(2)
matmjb(4)=eubeta3(beta1,eta,si)
matmjb(5:6)=(/0.0d0,0.0d0/)
end subroutine meanjbc

subroutine varjbc(matvjb,beta,eta,si)
!Variance of the moment conditions of JB
use ghs
use testmomsrot
use linear_operators
double precision matvjb(6,6),matmjb(6,1),beta(3),eta,si,beta1,nu,gam&
&,m400,m004,m044,m008,m800,m440
integer i,j

forall(i=1:6,j=1:6)
	matvjb(i,j)=0.0d0
end forall
beta1=dsqrt(dot_product(beta,beta))
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si

call meanjbc(matmjb(:,1),beta,eta,si)

m400=eubeta4(beta1,eta,si)
m004=ehr(nu,gam,2)*3.0d0
m044=ehr(nu,gam,4)*9.0d0
m008=ehr(nu,gam,4)*105.0d0
m800=eubeta8(beta1,eta,si)
m440=eub4hn(beta1,eta,si,2)*3.0d0    !E(e1^4*ej^4)

matvjb(1,1)=m800-6.0d0*m400+9.0d0
matvjb(1,2)=m440-3.0d0*m400-3.0d0*m004+9.0d0
matvjb(1,3)=m440-3.0d0*m400-3.0d0*m004+9.0d0
matvjb(2,1)=matvjb(1,2)
matvjb(2,2)=m008-6.0d0*m004+9.0d0
matvjb(2,3)=m044-6.0d0*m004+9.0d0
matvjb(3,1)=matvjb(1,3)
matvjb(3,2)=matvjb(2,3)
matvjb(3,3)=m008-6.0d0*m004+9.0d0

matvjb(1,4)=eubeta7(beta1,eta,si)-3.0d0*eubeta3(beta1,eta,si)
matvjb(2,4)=eub3hn(beta1,eta,si,2)*3.0d0-3.0d0*eubeta3(beta1,eta,si)
matvjb(3,4)=matvjb(2,4)

matvjb(4,4)=eubeta6(beta1,eta,si)
matvjb(5,5)=15.0d0*ehr(nu,gam,3)
matvjb(6,6)=matvjb(5,5)

forall(i=1:3,j=4:6)
	matvjb(j,i)=matvjb(i,j)
end forall

matvjb=matvjb-(matmjb.xt.matmjb)
end subroutine varjbc

subroutine var0u(wmat,dmut,dvsig,pth,sigma,b)
!Variance of the moment conditions under the null
use linear_operators
use matrixfuns
use datahadler
integer, parameter:: n=3,p=6
integer pth
double precision wmat(p,p),wmatmc(p,p),wmatu(pth+p,pth+p),dvsig(pth,n*n),dmut(pth,n)
double precision:: sigma(n,n),isigma(n,n),rsigma(n,n),irsigma(n,n),dmach,tol,b(n),beta(n)&
&,matV(n,n),irsigmah(n,n)&
&,matmsk(n,p),matvsk(n*n,p),matD(p,pth+p)
integer i,j,irank
!!!

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)

if (dot_product(beta(2:n),beta(2:n))<1.0d-10) then
	matV=eye(n)
else
!	Sigma=eye(n)
   matV(1,:)=(/-0.57735026918963d0, 0.81649658092773d0,0.0d0/)
   matv(2,:)=(/-0.57735026918963d0,-0.40824829046386d0,-0.70710678118655d0/)
   matv(3,:)=(/-0.57735026918963d0,-0.40824829046386d0, 0.70710678118655d0/)

!	matV(:,1)=beta/dsqrt(dot_product(beta,beta))
!	matV(1,2:3)=(/-0.50497514502648d0,-0.58533678384716d0/)
!	matV(2,2:3)=(/0.78236638609871d0,-0.0d0/)
!	matV(3,2:3)=(/-0.36455855607615d0,0.81079026232156d0/)

!	print*, "Please, supply orthogonal matrix!"
!	print*, "beta", beta
!	stop
	!matV(:,1)=beta/dsqrt(dot_product(beta,beta))
end if
irsigmah=matV.tx.irsigma

forall(i=1:p,j=1:p)
	wmatu(pth+i,pth+j)=0.0d0
end forall
forall(i=1:3)
	wmatu(pth+i,pth+i)=96.0d0
	wmatu(pth+3+i,pth+3+i)=15.0d0
end forall

forall(i=1:n,j=1:p)
	matmsk(i,j)=0.0d0
end forall
matmsk(1,4)=3.0d0
matmsk(2,5)=3.0d0
matmsk(3,6)=3.0d0
forall(i=1:n*n,j=1:p)
	matvsk(i,j)=0.0d0
end forall
matvsk(1,1)=12.0d0
matvsk(5,2)=12.0d0
matvsk(9,3)=12.0d0

wmatu(1:pth,pth+1:pth+p)=matmul(dmut.xt.irsigmah,matmsk)&
&+.5d0*matmul(dvsig.xt.kron(irsigmah,irsigmah),matvsk)

wmatu(1:pth,1:pth)=matmul(dmut,matmul(isigma,.t.dmut))&
&+.5d0*matmul(dvsig,matmul(kron(isigma,isigma),.t.dvsig))

forall(i=1:pth,j=pth+1:pth+p)
wmatu(j,i)=wmatu(i,j)
end forall

call exportvec(vec2(wmatu(pth+1:pth+p,1:pth)),"matB0.txt",2)

matD(:,pth+1:pth+p)=eye(p+1)
matD(:,1:pth)=-matmul(wmatu(pth+1:pth+p,1:pth),.i.wmatu(1:pth,1:pth))

wmat=(matD.x.wmatu).xt.matD
end subroutine var0u


end module moments3jb